#!/usr/bin/python
'''
Takes a SAM file of mapped assembled contigs and screens them for 50%
alignment against reference. Generates an AGP for it.
'''
import string
import sys
import fasta_read

# parameters.
sam_file = sys.argv[1]
contig_file = sys.argv[2]
agp_file = sys.argv[3]
fixed_file = sys.argv[4]

# read fasta file.
contigs = fasta_read.read_fasta_dict(contig_file)

# read sam file.
fin = open(sam_file, "rb")
lines = fin.readlines()
fin.close()

# parse hits.
written = {}
fout = open(agp_file, "wb")
fixed = open(fixed_file, "wb")
for line in lines:
	# tokenize.
	tmp = line.strip().split()
	cigar = tmp[5]
	
	
	# parse cigar.
	clist = []
	num = ""
	for c in cigar:
		# determine if number.
		try:
			i = int(c)
		except:
			i = None
	
		# not number.
		if i == None:
			# add pair.
			clist.append((c, int(num)))
			num = ""
			continue
			
		# is a number.
		num += c

	# get match length.
	mlen = 0
	for chit in clist:
		if chit[0] == "M":
			mlen += chit[1]
			
	# check if its greater than 50%.
	sz = len(contigs[tmp[0]])
	if float(mlen) > (float(sz) * .5):
		
		# check for duplicates.
		if tmp[0] in written:
			print "duplicate"
			print tmp
			sys.exit()
		
		# get orientation.
		if tmp[1] == "0":
			o = "+"
			fixed.write(">%s\n%s\n" % (tmp[0], contigs[tmp[0]]))
		else:
			o = "+"
			t = contigs[tmp[0]].upper()
			t = t.translate(string.maketrans("ATCGN", "TAGCN"))[::-1]
			fixed.write(">%s\n%s\n" % (tmp[0], t))
		
		# write agp.
		txt = "%s\t%i\t%i\t%i\t%s\t%s\t%i\t%i\t%s\n" % \
			("chr21", 1, sz, 1, "W", tmp[0], 1, sz, o)
		fout.write(txt)
fout.close()
fixed.close()
